{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "heading_collapsed": "false"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scanpy as sc\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.metrics import silhouette_samples, silhouette_score\n",
    "\n",
    "import sys\n",
    "sys.path.append('~/SCALEX/script')\n",
    "from metrics import *\n",
    "plt.rc('font', family='Helvetica')\n",
    "plt.rcParams['pdf.fonttype'] = 42"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "palette = { 'SCALEX': (0.8392156862745098, 0.15294117647058825, 0.1568627450980392),\n",
    "            'Raw': (0.09019607843137255, 0.7450980392156863, 0.8117647058823529),\n",
    "            'Seurat_v3': (1.0, 0.4980392156862745, 0.054901960784313725),\n",
    "            'Harmony':  (0.12156862745098039, 0.4666666666666667, 0.7058823529411765),\n",
    "            'Conos':(0.4980392156862745, 0.4980392156862745, 0.4980392156862745),\n",
    "            'BBKNN': (0.5803921568627451, 0.403921568627451, 0.7411764705882353),\n",
    "            'Scanorama': (0.8901960784313725, 0.4666666666666667, 0.7607843137254902),\n",
    "            'FastMNN':  (0.17254901960784313, 0.6274509803921569, 0.17254901960784313),\n",
    "            'scVI': (0.5490196078431373, 0.33725490196078434, 0.29411764705882354),\n",
    "            'online_iNMF':(0.7372549019607844, 0.7411764705882353, 0.13333333333333333),\n",
    "            'LIGER':(0.6509803921568628, 0.807843137254902, 0.8901960784313725)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets =['pancreas', 'heart','liver','NSCLC','PBMC']\n",
    "methods = ['Raw', 'SCALEX','Seurat_v3','Harmony','Conos','BBKNN','scVI','Scanorama','FastMNN','online_iNMF','LIGER']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": "true",
    "tags": []
   },
   "source": [
    "### overcorrection_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy\n",
    "from sklearn.neighbors import NearestNeighbors, KNeighborsRegressor\n",
    "def overcorrection_score(emb, celltype, n_neighbors=100, n_pools=100, n_samples_per_pool=100, seed=124):\n",
    "    \"\"\"\n",
    "    \"\"\"\n",
    "    n_neighbors = min(n_neighbors, len(emb) - 1)\n",
    "    nne = NearestNeighbors(n_neighbors=1 + n_neighbors, n_jobs=8)\n",
    "    nne.fit(emb)\n",
    "    kmatrix = nne.kneighbors_graph(emb) - scipy.sparse.identity(emb.shape[0])\n",
    "\n",
    "    score = 0\n",
    "    celltype_ = np.unique(celltype)\n",
    "    celltype_dict = celltype.value_counts().to_dict()\n",
    "    \n",
    "    N_celltype = len(celltype_)\n",
    "\n",
    "    for t in range(n_pools):\n",
    "        indices = np.random.choice(np.arange(emb.shape[0]), size=n_samples_per_pool, replace=False)\n",
    "        score += np.mean([np.mean(celltype[kmatrix[i].nonzero()[1]][:min(celltype_dict[celltype[i]], n_neighbors)] == celltype[i]) for i in indices])\n",
    "\n",
    "    return 1-score / float(n_pools)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oc_score = pd.DataFrame(index = datasets, columns = methods)\n",
    "\n",
    "for dataset in datasets:\n",
    "    adata = sc.read_h5ad('~/SCALEX/results/{}/adata.h5ad'.format(dataset))\n",
    "\n",
    "    for method in methods:\n",
    "        if method+'_umap' in list(adata.obsm.keys()):\n",
    "            oc_score.loc[dataset, method] = overcorrection_score(adata.obsm[method+'_umap'], adata.obs['celltype'])\n",
    "        else:\n",
    "            oc_score.loc[dataset, method] = 0\n",
    "            pass\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5116bdb5e0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "findfont: Font family ['Helvetica'] not found. Falling back to DejaVu Sans.\n",
      "findfont: Font family ['Helvetica'] not found. Falling back to DejaVu Sans.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoMAAAGDCAYAAABQhhoTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABCQUlEQVR4nO3dd7icRdnH8e+PLgkhlIj0IAKhI4KIVAFRmiC+0ktUREAEX7qAckBEQEBepIMSaqR3kR6p0qSG0AIJCSESSiqhJNzvHzMLD8vuOXtyzu4p+/tc117ZnZln9t7lcHJn5pkZRQRmZmZm1pxm6+oAzMzMzKzrOBk0MzMza2JOBs3MzMyamJNBMzMzsybmZNDMzMysiTkZNDMzM2tiTgbNzMzMmli3SAYlLSjpeknTJI2WtEuVdodKek7SFEmvSTq0rH6UpOmSpubHHY35BGZmZmY90xxdHUB2FvARsAiwBnCrpKcjYnhZOwF7AM8AywJ3SBoTEX8vtNkmIu5qQMxmZmZmPZ66+gQSSX2A94BVIuKlXHYp8EZEHNHGtWeQPsOv8utRwF7tTQYXXnjhGDhw4CxEb2Zm1vM88cQTb0fEgK6Ow7qH7jAyuDwwo5QIZk8DG7V2kSQBGwDnlVVdLmk24Eng0Ih4usr1ewN7Ayy11FI8/vjjsxi+mZlZzyJpdFfHYN1Hd7hnsC8wuaxsEjBfG9e1kOK/qFC2KzAQWBq4F7hdUv9KF0fE+RGxVkSsNWCA/3FkZmZmzak7JINTgX5lZf2AKdUukLQ/6d7BrSLiw1J5RDwYEdMj4v2I+CMwkTR6aGZmZmYVdIdk8CVgDknLFcpWB8oXjwAg6afAEcCmETG2jb6DtOjEzMzMzCro8mQwIqYB1wHHSeojaT1gW+DS8raSdgVOAL4bEa+W1S0laT1Jc0maJ287szDwYP0/hZmZmVnP1OXJYLYf8CXgLWAosG9EDJe0gaSphXbHAwsBjxX2Ejw3180HnENamfwG8H1gi4h4p2GfwszMzKyH6Q6riYmId4HtKpTfT1pgUnq9TCt9DAdWq0d8ZmZmZr1VdxkZNDMzM7Mu4GTQzMzMrIk5GTQzMzNrYk4GzczMzJqYk0EzMzOzJuZk0MzMzKyJORk0MzMza2LdYp9BMzPrPe6+Z9mqdZtuMrKBkZhZLTwyaGZmZtbEnAyamZmZNTEng2ZmZmZNzPcMmlnT8L1sZmZf5JFBMzMzsybmZNDMzMysiTkZNDMzM2tiTgbNzMzMmpiTQTMzM7Mm5mTQzMzMrIk5GTQzMzNrYk4GzczMzJqYk0EzMzOzJuZk0MzMzKyJORk0MzMza2JOBs3MzMyamJNBMzMzsybmZNDMzMysiTkZNDMzM2tiTgbNzMzMmpiTQTMzM7Mm5mTQzMzMrIk5GTQzMzNrYk4GzczMzJrYHF0dgJlZZxr65NiqdV9uYBxmZj2FRwbNzMzMmpiTQTMzM7Mm5mTQzMzMrIk5GTQzMzNrYl5AYmZm7dLaIh3wQh2znsYjg2ZmZmZNzMmgmZmZWRNzMmhmZmbWxJwMmpmZmTUxJ4NmZmZmTczJoJmZmVkTczJoZmZm1sScDJqZmZk1MSeDZmZmZk3MyaCZmZlZE3MyaGZmZtbEnAyamZmZNTEng2ZmZmZNzMmgmZmZWRNzMmhmZmbWxJwMmpmZmTUxJ4NmZmZmTaxbJIOSFpR0vaRpkkZL2qVKu0MlPSdpiqTXJB1aVj9Q0r2S3pf0gqTNGvMJzMzMzHqmbpEMAmcBHwGLALsC50hauUI7AXsACwDfB/aXtFOhfijwJLAQcBRwjaQB9QzczMzMrCfr8mRQUh/gR8BvI2JqRDwA3ATsXt42Ik6OiP9ExIyIeBG4EVgv97M8sCZwTERMj4hrgWdz32ZmZmZWQZcng8DywIyIeKlQ9jRQaWTwU5IEbAAMz0UrA69GxJRa+pG0t6THJT0+YcKEWQ7ezMzMrCfrDslgX2ByWdkkYL42rmshxX9RoZ9JtfYTEedHxFoRsdaAAZ5JNjMzs+Y0R1cHAEwF+pWV9QOmVGgLgKT9SfcObhARH85qP2ZmZmbNrjuMDL4EzCFpuULZ6nw2/fs5kn4KHAFsGhFjC1XDga9KKo4EVu3HzMzMzLrByGBETJN0HXCcpL2ANYBtgW+Xt5W0K3AC8J2IeLWsn5ckPQUcI+loYAtgNbyAxMx6iIFH3Fq1btSJWzUwEjNrJl2eDGb7AX8D3gLeAfaNiOGSNgBui4i+ud3xpG1jHkvrRwC4LCL2yc93AoYA7wGvA/8TEV4dUqOhT46tWrfz15doYCRmZmbWKN0iGYyId4HtKpTfT1oYUnq9TBv9jAI27tzozMzMzHqv7nDPoJmZmZl1ESeDZmZmZk3MyaCZmZlZE3MyaGZmZtbEnAyamZmZNTEng2ZmZmZNzMmgmZmZWRNzMmhmZmbWxJwMmpmZmTUxJ4NmZmZmTczJoJmZmVkTczJoZmZm1sScDJqZmZk1sTna01jSwsCywFMR8WF9QrLu6O57lq1at+kmIxsYiZmZmXWmmkYGJc0n6SrgLeAhYPFcfq6klvqFZ2ZmZmb1VOs08UmkBHBNYHqh/Bbgh50dlJmZmZk1Rq3TxD8AfhgRT0mKQvkI4KudH5aZmZmZNUKtI4MLAO9UKJ8PmNl54ZiZmZlZI9WaDD5GGh0sKY0O/oJ0D6GZmZmZ9UC1ThMfCdwuaeV8zUH5+TeBDesVnJmZmZnVV00jgxHxELAuMBcwEtgUGAesGxH/qV94ZmZmZlZPbY4MSpoTuAw4MiL2rH9IZmZmZtYobY4MRsTHwOZ8dp+gmZmZmfUStS4guQ7Yvp6BmJmZmVnj1bqA5HXgaEkbAI8D04qVEXFaZwdmZmZmZvVXazI4GHgPWC0/igJwMmhmZmbWA9WUDEbEMvUOxMzMzMwar9Z7Bj8lqa+kPvUIxszMzMwaq9ZpYiT9EjgcWDy/HgucFBFn1yk2MzPLVr141ap1z+75bAMjMbPepqZkUNKRwG+AU4AHcvEGwImS+kXEiXWKz8zMzMzqqNaRwX2AvSNiaKHsbkkvAycATgbNzMzMeqBa7xn8MvBYhfJHgUU6LxwzMzMza6Rak8GXgF0qlO8CvNh54ZiZmZlZI9U6TdwCXCVpQ+DBXLYesBHw4zrEZWZmZmYNUNPIYERcB6wDjAe2zo/xwDcj4oa6RWdmZmZmdVXz1jIR8QSwWx1jMTMzM7MGq2lkUNKPJW1boXxbSf/T+WGZmZmZWSPUuoCkBfigQvm0XGdmZmZmPVCtyeBXqbxq+JVcZ2ZmZmY9UK3J4HvAchXKlwemdF44ZmZmZtZItSaDNwJ/lrR8qUDSCsBpwA11iMvMzMzMGqDWZPBwYBLwvKQxksYAw4HJwKH1Cs7MzMzM6qumrWUiYjKwnqTvAmvk4ieBuyMi6hSbmZmZmdVZzfsMAkTEncCddYrFzMzMzBqs1WliSatL+k5Z2a6SXpX0lqRzJc1V3xDNzMzMrF7aumfweGD90gtJKwEXAS8DQ4FdSfcTmpmZmVkP1FYyuCZwR+H1TsDzEfG9iDgQ+DWwY51iMzMzM7M6aysZXAgYV3i9IXBz4fUwYKlOjsnMzMzMGqStBSQTgMWBMZJmB74BnFKonwv4pE6xmZmZWTczYtCKdd1FZMUXRqie/dsXtTUyOAw4RtJXgYNz2b2F+pWAUZ0flpmZmVnbJK0v6SFJkyS9K+lBSWvnukUl/VXSm5KmSHpB0rGS+hSuV14Y+3yFvodJ2qtC+UBJIWlq2WPHXH+tpAvKrrle0pmd/w10XFsjg78F7iKdQTwTOCAiphXqdwfurlNsZmZmZlVJ6gfcAuwLXEWasdwA+FDSgsDDwEPAuhExStKSwCHAssAzuZsNgS8Dc0haOyIea0cI/SNiRoXyXwLDJV0REffmJHFNUt7U7bSaDOYvbhCwMjAhIsaVNTkGGFuv4MzMzMxasTxARAzNr6eTF75KOh6YAuwWEZ/kdmOAA8v62JN07O6X8vP2JIMVRcR4SQcDF0jaBDgD2DMipna073po8zi6iJgREU9XSATJ5e/UJzQzMzOzVr0EzJR0saQtJC1QqNsMuK6UCFYiaV7gf4DL82Onzto/OSKGACOB/wD/jIh/dka/9VDr2cRmZmZm3Uo+Lnd9IIALgAmSbpK0CGlHlDfb6GJ74EPSaOKtwJzAVu0I4W1JEwuPFcvq789xXNaOPhvOyaCZmZn1WBExIiIGR8QSwCrAYsDpwDvAom1cvidwVZ4F/QC4NpfVauGI6F94jChVSFqOdH/i2cCpkuZsR78N1S2SQUkL5lU20ySNlrRLlXbfkXRvXjE0qkL9KEnTC6t67qjQjZmZmfVCEfECMISUFN4F/FBSxVxH0hLAJsBuksZLGk+aMt5S0sIdiUOSgAtJSemvgGl04xPbukUyCJwFfAQsQjri7hxJK1doNw34G3BoK31tExF982Pzzg/VzMzMugNJgyQdnBM78mrhnYF/A6cB/YCLJS2d6xeXdJqk1Ugre18CVgDWyI/lSQtjdy68zRyS5ik8ahnh2xdYGDgh37P4M+CwvCi322lraxkAJM0EFo2It8rKFwLeiojZZzWAvNfPj4BV8iqbByTdRPqPdESxbUQ8CjwqabNZfT+zrjT0yeqL73f++hINjMTMbNZ0s02hpwDrAAdJ6g9MJG01c2hETJb0beB44JGcb7wBDCVtmbcncFZEjC92KOncXPeXXHROfpRcDhydn09Mg4Cf+h1wDXAC8P2I+AggIp6XdCppdfGGEVHXjbvbq6ZkEKj2H35u0oheRywPzIiIlwplTwMbzWJ/l+ch4SdJPwxPV2okaW9gb4CllvKJemZmZj1NRLwB7NBK/Tjgp1WqK47SRcTJwMn5+catvH1rSXH/Cv0eCxzbyjVdptVkUNJB+WkA+0gq7o8zO2ljxxc6GENfYHJZ2SRgvlnoa1fSEm6R9hG6XdKgiJhY3jAizgfOB1hrrbW6VYZuZmZm1ihtjQz+Kv8pYC/SKSQlH5GOotungzFMJc3pF/UjDf22S0Q8WHj5R0l7khLWm2c9PDMzM7Peq60TSJYBkHQvsH1EvFeHGF4i3Zy5XES8nMtWB4Z3Qt9B68O4ZmZmZk2tptXEEfGdOiWC5LOOrwOOk9RH0nrAtsCl5W0lzSZpHtKmkMqreubKdUtJWk/SXLn8UNJKngfL+zEzMzOzpNYFJORDljclHeb8uSQyIn7QwTj2I20Z8xZpk8h9I2K4pA2A2yKib263IXBv4brpwL+AjUn3GJ5DOnz6A+ApYAsfl2dmZmZWXa1by/wJ+DUpERtHmn7tNBHxLrBdhfL7SQtMSq+HUWXaNyKGA6t1ZlxmZmZmvV2tI4N7ADtHxDX1DMbMzMzMGqvWE0hmI027mpmZmVkvUuvI4PnAbkBL/UIxMzOzbq9l/vruzdsyybuANFityWB/YBdJ3wWeAT4uVkbEAZ0cl5k1Wsv8rdRNalwc1m4tLS2zVGfWG0han3RiyMqk/ZBHAL+OiMcaHMcoYK+IuKuD/dwLrEI65e014HcRcWPHI6yu1mRwJT6bJi4/vsWnd5iZmVnDSepHOot4X+AqYC7SYRMfdvL7zBERMzqzz1YcCDwfETMkrQPcJWn5iHizXm/Ynn0Gqz02qVdwZmZmZq1YHiAihkbEzIiYHhF3RMQzAJJ+KmmEpPck3S5p6Vw+UFJI+nRQTNIwSXvl54MlPSjpz5LeAVokLSvpHknvSHpb0uWS+uf2lwJLATdLmirpsGoBS7pN0v5lZU9L2j5/lmcKiWeQ9lZeslO+rSpqXUACQN7MeRVJK+fNn83MzMy6ykvATEkXS9pC0gKlCknbAkcC2wMDgPuBoe3oex3gVWAR4A+kre3+CCwGrEhK0FoAImJ34HVgm4joGxEnt9LvUGDnQpwrAUsDtxbKbpH0AfAIMAx4vB1xt1tNyaCkOfNeg+8BTwPPAu9JOlnSnPUM0MzMzKySiJgMrE8aQbsAmCDpJkmLAPsAf4yIEXmk7QRgjdLoYA3GRcRfImJGHnF8JSLujIgPI2ICcBqw0SyEfX1ZHLsC10XEp1PbEbE16TCNLYE7IuKTWXifmtU6MngSaTXxPqQh2eVI8/O7k7JkMzMzs4bLyd7giFiCtPBiMeB00mjb/0maKGki8C5pdG/xGrseU3whaRFJf5f0hqTJwGWkY2/bG+8U0ijgTrloZ+DyCu0+jojbgM0ldfSkt1bVmgzuAvwsIi6OiJH5MQTYi5TRmpmZmXWpiHgBGEJKCscAv4iI/oXHlyLiIWBavmTewuVfKe+u7PUJuWzViOhHGiRTK+1bMxTYWdK6wDx8/qjdcnOQjtqtm1qTwfmBkRXKR5K2nTEzMzNrKEmDJB0saYn8eknSSNu/gXOB30haOdfNL+nHAHma9w1gN0mzS/opbSdc8wFTgUmSFgcOLav/L/DVGkP/B2nk8jjgytI0cP48W0j6Ur5FbzdgQ+BfNfY7S2rdWuZp4ADgl2XlB+KTSczMzJpH99oUegppocdBeWXvRNJWM4dGxGRJfYG/5/vzJgF3Alfna38OnE0a8fsr8FAb73UscEnu5xXgUuB/C/V/BP4i6WTg+Ig4pVpHEfGhpOuAn5IWuZSItChlJdKeiS8DO0bEf9qIrUNqTQYPA/4haTNStg3wLdK8/Bb1CMzMzMysNRHxBrBDK/WXkpK2SnW3ActUqRtCmm4ulg0HvlHW9NRC/Y1AzZtDR8TPgJ+VlY0gJbcNVes+g/eRFo5cA/TNj6uBFSLigfqFZ2ZmZmb1VOvIIBExDjiqjrGYmZmZ9XiSdgXOq1A1OiJWbnQ8bamaDEpaE3gqIj7Jz6uq91y2mZmZWU8REZdTYbuY7qq1kcHHScus38rPg88voS4JYPbOD83MzMzM6q21ZHAZYELhuZnV0d33VN/VYNNNKu3sZGZm1nFVk8GIGF18CYyJiC9sqChpqXoEZmZmZmb1V+um06+RDnn+HEkL5TozMzMz64FqTQZF5WNW+gIfdF44ZmZmZtZIrW4tI+mM/DSAP0p6v1A9O/BNfAKJmZlZ0xh4xK3tOYO33UaduFV3OuGkKbQ1MrhqfghYsfB6VeBrwH+AwXWMz8zMzKwiSaPy6WjFssGSfCBGO7Q6MhgR3wGQdBFwYERMbkhUZmZmZg0kaY6ImNHVcXSFWu8Z/A3Qr7xQ0hKSFunckMzMzMw6TtIRkkZKmiLpeUk/LNQNlvSgpD9LegdokTRE0tmSbpM0Ndd/RdLpkt6T9IKkrxf6WFHSMEkTJQ2X9INC3RBJZ0m6Nb//I5KWzXVnSTq1LNabJP1vA76WL6g1GbwM2KJC+feocgC0mZmZWRcbCWwAzA8cC1wmadFC/TrAq8AiwB9y2Q7A0cDCwIfAw6Tb4hYGrgFOA5A0J3AzcAfwZeBXwOWSVij0v1N+3wWAVwrvcTGws6TZcl8LA5sBV3TS526XWpPBtYD7KpTfn+vMzMzMusINeWRuoqSJwNmlioi4OiLGRcQnEXEl8DJp8WvJuIj4S0TMiIjpuez6iHgiIj4Argc+iIhLImImcCVQGhn8FmlXlRMj4qOIuAe4Bdi50P/1EfFonn6+HFgjx/UoMAnYNLfbCRgWEf/tpO+kXWpNBucA5q5QPk+VcjMzM7NG2C4i+pcewH6lCkl7SHqqkCiuQhrhKxlTob9iQja9wuu++flipAM5PinUjwYWL7weX3j+fuFaSKODu+Xnu9GFM621JoOPAPtWKP8l8FjnhWNmZmbWcZKWBi4A9gcWyonic6QdUko6sk3OOGDJ0lRvthTwRo3XXwZsK2l10o4tN3Qglg5pdTVxwVHAPZJWA+7JZZuQhko3q3qVmZmZWdfoQ0r2JgBI+glpZLCzPEIa7TssLwZZD9gGWLuWiyNirKTHSCOC1xamqRuupmQwIv4taV3gMGD7XPwksF9EPF2v4MzMepWW+dto0CX3jpu1S0/ZFDoins9J2sPAJ8AlwIOd2P9HkrYh3aP4G9KI4B4R8UI7urmYlAwe2FlxzYpaRwbJSd+udYzFzMzMrGYRMbBC2RBgSH5+FGl2s9K1n7YrlA0ue30hcGHh9SsUcqeIGA5sVKX/8r6GAUuUNXuddN/isEp9NEqt9wwiaRFJh+T9dxbOZetJWqZ+4ZmZmZn1PnlrmgOBCyOirkf8taWmZFDSN4AXSSODe/HZBtTf5bM9c8zMzMysDZJWBCYCiwKnd2kw1D4yeArwfxHxddIGjCW3k26YNDMzM7MaRMSIiOgTEd/uDkf91poMfoN0k2O5N0m7dpuZmZlZD1RrMjiddJRKuUHAW50XjpmZmZk1Uq3J4I3AMZJKp42EpIHAScC19QjMzMzMzOqv1mTwEGBB0saN8wIPkA5cnkg6zNnMzMzMeqBa9xmcAWwMbAisSUoi/xMRd9UpLjMzMzNrgDaTQUmzA5OA1SPiHj47js7MzMyazKoXr1rXPfGe3fPZHnHCSW/SZjIYETMljQbmakA8ZmbWic7ap/q/33957iYNjMSs80kaRdrVZGahePmIGNfOfgYDe0XE+oWyIcCewHYRcWOh/M/Ar4GfRMSQfO1FwOERcXKh3Vhgt4gYJqkFOAbYMSKuyvVzAB8Dy0TEqPbE29lqvWfw98CJpZNHzMzMzLqJbSKib+HRrkSwDS8Be5Re5ARuB2BkWbt3gcMkzddKX+8Cx+YZ126lPQtI1gfekDRS0jPFRx3jMzMzM6uJpAUk3SJpgqT38vMlCvWDJb0qaYqk1yTtmk8DORdYV9JUSRMLXd4MrC+ptL3e94FngPFlbz0CeBg4qJXw/gl8BOzWsU/Z+WpdQHJNXaMws16rpaVllurMzGbBbKQp2x2A2YG/AWcC20nqA5wBrB0RL0paFFgwIkZI2oeyaeLsA9L2ejsB55BGCS8BflnhvX8L3CvpLxHxboX6yG1Ol3RFft0t1LKAZE6gD3BWRIyuf0hmZmZmNbtB0oz8fFhEbFeqkPQH4N5C20+AVSS9HhFvkk5Sa8slwJ8kDQU2It1H+IVkMCKeknQncHh+fEFE3CTpKGAv4IIa3rsh2pwmjoiPgX0Br+4xMzOz7ma7iOgfEf2BXSSdJ2m0pMnAfUB/SbNHxDRgR2Af4E1Jt0oa1FbnEfEAMAA4CrglIqa30vx3wL6SWjuq9+jc1zw1fboGqPWewTsALzszMzOz7uxgYAVgnYjoR9ofGfKAVkTcHhHfBRYFXuCz0bm2pmwvy31f0lqjiHgBuI6U7FVrcyfp4I792njPhqn1nsG7gRMkrQY8AUwrVkbEdZ0dWG8x8Ihbq9aNOnGrBkZiZmbW680HTAcmSlqQtJ0LAHm07lvAXbnNVNK0McB/gSUkzRURH1Xo9wzgftJIY1uOJS0yaW1G9SjSvYjdQq3J4Jn5zwMq1AXpJk0zMzPr5br5ptCnA1cAbwPjgFOB7XLdbKTVvpeQcpenSLfBQTpQYzgwXtInEfG5rfTygpC7awkgIl6TdGmh70ptHpT0KLBFLX3WW03JYETUOp1sZmZm1hARMbDs9TjS8blF5+U/3yQtAKnUz0fAVmVlg1t53/ULz4cAQ8rq96MwDRwRLRX62LJa/43mJM/MzMysidU6TYykrUhLpVciDa8+D5wUEf+oU2xWhY+XMjMzs85SUzIoaS/gbOBy4OJcvAFwvaR9I+JvdYrPzLoBL4QyM+u9ah0ZPBw4KCLOLJT9VdITwBGkHb7NzMzMrIep9Z7BpUhn6pW7DVi6o0FIWlDS9ZKm5Y0id6nS7juS7pU0SdKoCvUDc/37kl6QtFlHYzMzMzPrzWpNBl8HvluhfHOgM46oO4t0ePMiwK7AOZJWrtBuGmkU8tAq/QwFngQWIu3hc42kAZ0Qn5mZmVmvVOs08SnAXyStCTyUy9YDdgd+1ZEA8sHRPwJWiYipwAOSbsp9H1FsGxGPAo9WGvGTtDywJrB5PirmWkm/zn2f25EYzczMzHqrWvcZPE/SW6SjWLbPxSOAHSKioztoLw/MiIiXCmVPU2UvoFasDLwaEVPK+qk0woikvYG9AZZaaql2vpWZmZlZ71Dz1jIRcT1wfR1i6AtMLiubRDpSpr39TKrQz+KVGkfE+cD5AGuttVZbZxKamfU6Y4+4v2rdla+dVLVusSM82dLMWlpa6vp3ZktLS3c+4aRXqumeQUkbSfrCSF0u37DSNe0wFehXVtYPmFKhbSP6MTMzsx5C0i6SHpc0VdKbkm6TtH7bV1pJrQtI/gwsUKG8X67riJeAOSQtVyhbnXRGYHsMB74qqTiiOCv9mJmZWQ8g6SDSecQnkBahLkXaF3nbLgyrx6l1mngF0v135Z7LdbMsIqZJug44Lm9uvQbpP+K3y9tKmg2YC5gzvdQ8wCcR8VFEvCTpKeAYSUeTDn9ejbSApMdpaWmpWjeAjg7GmpmZ9WyS5geOA34SEdcVqm4GbpY0N3ASsEMuvwo4PCI+lLQxcBlpQOtwYCZwZERcVOj7L6Rc4n3gAuCEiPhE0teAv5LylY+BuyNixzp+1LqrNRmcDiwKvFZWvjhpS5iO2o+0ZcxbwDvAvhExXNIGwG0R0Te32xC4tyyuf/HZodQ7kQ6Lfo+0Hc7/RMSETojP7HNau9dqiRM3aGAkZmZNa11gHqqvZzgK+BYpaQvgRuBo4Le5/ivA/KRc5ruk7ehuiIj3SIng/MBXSdvV3QG8SUoCf59ff4c0QLVWJ3+uhqs1GbwdOEnSD/KXhKQFgT/mug6JiHeB7SqU309aGFJ6PQyoemNpRIzis8SwcVrmb6XyioaFYWZm1kQWAt6OiBlV6ncFfhURbwFIOhY4j8+SwY+B4/L1/5A0FVhB0mOkwaU18g4lUySdStry7q/5uqWBxSJiLPBAfT5e49SaDB4C3AeMkvRMLluNNJLXo4dGu9KqF69ate5HPXN228zMrFHeARaWNEeVhHAxPn8wxuhc9un1Zde9TxqAWph0O1r5taXdSQ4jjQ4+Kuk94NSI6NHH8ta0gCQi3iQtxjgEeCY/DgZWj4hx9QvPzMzMrKKHgQ+pMLOYjePzR+Yulcva8jafjf4Vr30DICLGR8TPI2Ix4BfA2fk+wh6rPfsMlm6gNDPrUqfuuHXVOu+BZ9YcImKSpN8BZ0maQbqP72NgM9L9fEOBo/O0bwC/Iy0aaavfmZKuAv4gaQ9gQeAg0mlsSPox8HCeIn4v9/1JZ3++Rqo5GSyRNJk0j/5qHeKxOmrtL1DwX6JmZta27rQpdEScKmk8aWHI5aS9hZ8A/gD8h7QFXun2tquB42vs+lekRSSvAh+QBsNKU8FrA6fnFcf/BQ7s6TlRu5NBWlnAYWZmZtZIEXE5KRGs5ID8KL9mGLBEWdnAwvP3gN2qvN9hpPsGe41aN502MzMzs16ozWRQ0pySrpS0bC66jC+eJWxmZmZmPVCb08QR8bGkzYHf5Nf71j0qM2sKZ+1zT9W6X567SQMjMTNrXrVOE18HbF/PQMzMzMys8WpdQPI6aXn2BsDjwLRiZUSc1tmBmVnP0Nrm6eAN1M3Murtak8HBpL10VsuPogCcDFq31Fqi8uyezzYwEjMzs+6ppmQwIpapdyBmjdbS0jJLdWZmZr1Ju7eWkbSIJG9JY2ZmZtYL1DQyKGlO0m7e+wJfApYHXpV0EjA6Is6uX4hmZmbWXZy1zz1Rz/5/ee4mPtyiwWq9Z/AYYBvSbtxXFMofBQ4HnAx2E2OPuL+rQzAzM2sISaOARYCZpHOJHwL2iYgxkoYAuwAfkdY3vAQcFBH/ytcOBvaKiPXz637AbcB4YGfgfGBPYJ2IeDS3+RrwckQovx4GfAtYLiLG5LLNgAuLJ5p0d7VO9+5M+nJv5POHMT9HGiU0MzMz6wrbRERfYFHSWcF/KdSdnOv6AecA10mavbwDSQsAdwOjgR0j4qNc9S5tn2c8Dfhtxz5C16p1ZHAx0hdU6fpZOd/YrNc6dcetq9YtdsS5DYzEzKx5RMQHkq4BTq9QF5KuAC4gjSSOK9VJGgDcCTwJ/CwiioNeFwO7SNqoNKJYwRnAIZJOioiRnfNpGqvWkcHhwIYVyncAnui8cMzMzMzaT9K8wI7AvyvUzQ7sAbxGGj0sWRAYBjwM/LQsEQR4HziBtG6imjdISeaxsxp7V6t1VO9Y4DJJSwKzAz+WNIg0F79VvYIzMzMza8MNkmYAfYAJwPcKdYdI2h+YGxBp5G9moX5JYB5SIlhtYcx5uZ8tgJertPkj8IqklTvwObpMTSODEXEzaRRwc9I9g8cAy5Hm6e+qX3hmZmZmrdouIvqTkrr9gX9J+kquOyXXzQusBfwpJ3UlTwOHALdJ+nqlziPiQ+D3+VFRREwAzgSO69hH6Ro17xcYEbdHxEYR0Tci5o2I9SPijnoGZ2ZmZlaLiJgZEdeRVhavX1YXEfEc8CBlM5oR8X/AicCdklap0v1FQH9g+1ZC+BPwHeAbs/QBulBNyaCkGyT9SNJc9Q7IzMzMrL2UbAssAIyoUD+IlCQOL6+LiJOB/wPukrRChfoZpFnRw6u9f0RMBE4FDpvFj9Blar1n8H3SipqPJV0LXNrKqhozMzPrpbrhptA3S5pJ2ktwNLBnRAyXBHCYpF+T7hd8hzTCd16lTiLi95LmBu6WtFGFJkOB35AWnVTzf8CBs/pBukqtZxPvIqkP8EPSopE7Jb1J+mIuy0OvZmZmZg3T2sbOETEYGNxK/RBgSFnZ0cDR+eXgsrpPgFXKyjYuez0V+HKrQXdD7blncFpEXBYRWwKLk+bGtwaeqlNsZmZmZlZn7d4wWtI8wCakpdvLA2M6OygzMx+taGbWGLUuIJGkzSVdTNqs8RzS7t2bRsQy9QzQzMzMzOqn1pHBN0nn+t1GmkO/tXBun5mZmZn1ULUmg78Frs7Lps3MzMysl6h1NfEFkuaXtFYuesWJoZmZmVnP1+Y9g5KWknQzaX+eR/LjbUk3SVq63gGamZmZWf20OjIoaXHg36TziH8HPJ+rVgb2Ax6StHZEjKtrlN3AiEErVq1bcacGBmJmZmbWidqaJj4GeA3YLCKmF8pvkPRn4I7c5hd1is+sS5y1zz2t1m/bf84GRWJm1r2MPeL+qGf/S5y4QXc74aTXa2uaeEvgyLJEEICIeJ+0S/dWX7jKzMzMrAEkrS/pIUmTJL0r6UFJa3d1XD1JWyODA4CRrdS/ktuYmZmZNZSkfsAtwL7AVcBcwAbAh10ZV4mkOSJiRlfH0Za2RgbfAr7WSv1yuY2ZmZlZoy0PEBFDI2JmREyPiDsi4hkAST+XNELSFEnPS1ozlx8haWSh/IelDiUNlvSApFMkvSfpNUlbFOoXy4to35X0iqSfF+paJF0j6TJJk4HBkr4p6WFJEyW9KelMSXMVrglJ+0l6Ocfze0nL5tHOyZKuKrWXtICkWyRNyLHdImmJjn6JbSWDtwHHS5q7vCIfS/d74B8dDcLMzMxsFrwEzJR0saQtJC1QqpD0Y6AF2IN0cMYPSDujQJr13ACYHzgWuEzSooV+1wFeBBYGTgb+Kql0L+PfgbHAYsD/ACdI2qRw7bbANUB/4HJgJvC/ua91gU1Ji3CLvgd8A/gWcBhwPrAbsCSwCrBzbjcbcBGwNLAUMB04s5YvqjVtJYMtwFeBVyQdLmnb/PgN8DKwLHBcR4MwMzMza6+ImAysDwRwATAhj9otAuwFnBwRj0XySkSMztddHRHjIuKTiLiSlNN8s9D16Ii4ICJmAhcDiwKLSFoSWA84PCI+iIingAtJCWfJwxFxQ+57ekQ8ERH/jogZETEKOA/YqOyjnBwRkyNiOPAccEdEvBoRk0gDc1/Pcb8TEddGxPsRMQX4Q4W+2q3VewYjYpykbwNnAycApaw4gNuB/SPijY4GYdamlvlbqZvUuDjMzKxbiYgRpKNykTQIuAw4nTSqVnHdg6Q9gIOAgbmoL2nkrmR8of/386BgX2Ah4N2ciJWMBtYqvB5T9l7LA6flNvOScq8nykL6b+H59Aqvv5L7mhf4M/B9oDQKOp+k2XPiOkva3HQ6IkZFxJakL+lb+TEgIraMiFdn9Y3NzMzMOlNEvAAMIU2tjiHNYH5OPjDjAmB/YKGI6E8ajatlS5txwIKS5iuULQUUB8bKt945B3gBWC4i+gFH1vhelRwMrACsk/vaMJd3aDueNpPBkoh4LyIezY93O/KmZmZmZh0laZCkg0uLKPI07s6kAzMuBA6R9A0lX8uJYB9SwjYhX/MTUvLYpogYAzwE/FHSPJJWA35GGo2sZj5gMjA1j1zuOyuftdDXdGCipAVJez13WE1nE5t1ZwOPuLVq3XzVD44xM7NZ0M02hZ5CWuxxkKT+wETSVjOHRsRkSQsBVwCLA6OA3SPiSUmnAg+TTli7BHiwHe+5M3AuaZTwPeCYiLirlfaHkBaEHAY8CVwJbNJK+9acTvo8b+f3PxXYbhb7+pSTQTMzM+uR8rqFHVqpP5eUuJWXHwUcVeWaIaSp5mKZCs/HAltXubalQtl9wKCy4t9V6ju/Xr/s9dGF5+OAjcv6Oq9SLO1R8zSxmZmZmfU+TgbNzMzMmpiTQTMzM7Mm5mTQzMzMrIk5GTQzMzNrYk4GzczMzJqYk0EzMzOzJuZk0MzMzKyJedNpMzMzq9mpO25dfvZupzr4ylu60wknTcEjg2ZmZtYUJL0g6acVyg+U9Hh+PkzSXo2Prus4GTQzM7NmcTGwR4Xy3XNdU3IyaGZmZj2WpMMlvSFpiqQXJW0qaXZJR0oamcufkLQkcCmwvqSlC9evBKwGDO2qz9DVukUyKGlBSddLmiZptKRdqrSTpJMkvZMfJ0lSoT5yH1Pz48LGfQozMzNrJEkrAPsDa0fEfMD3gFHAQcDOwJZAP+CnwPsRMRa4lzQSWLI78I+IeLuBoXcr3WUByVnAR8AiwBrArZKejojhZe32BrYDVgcCuBN4DTi30Gb1iHil3gGbmZlZl5sJzA2sJGlCRIwCyPf8HRYRL+Z2TxeuuRj4LXC8pNmAXYEDGxdy99PlI4OS+gA/An4bEVMj4gHgJj6ftZfsCZwaEWMj4g3gVGBww4I1MzOzbiMP/vwaaAHekvR3SYsBSwIjq1x2HbCopG8BGwPzArfWPdhurMuTQWB5YEZEvFQoexpYuULblfl8dl+p3X2Sxku6TtLAam8qaW9Jj0t6fMKECbMYupmZmXWliLgiItYHlibNGp4EjAGWrdL+feAa0kKS3YG/R8RHDQq3W+oOyWBfYHJZ2SRgviptJ5W161u4b3AjYCAwCBgH3CKp4lR4RJwfEWtFxFoDBgzoQPhmZmbWFSStIGkTSXMDHwDTgU+AC4HfS1ourzdYTdJChUsvBnYkzUw27Sriku5wz+BU0s2dRf2AKTW07QdMjYgAiIj7cvlHkg4kJZkrAs92asRmZmZNqpttCj03cCLp7/qPgYdI6wv+m+vuABYGXgB+WLjuPtKA0gcR8VgjA+6OukMy+BIwh6TlIuLlXLY6UL54hFy2OvBoG+1KAuhOP7RmZmbWSSLiGeCbVaqPz49K1wXw1Sp1G3dKcD1Il08TR8Q00s2cx0nqI2k9YFvSXkDlLgEOkrR4vkH0YGAIgKSVJa2R9xbqS1pc8gYwohGfw8zMzKwn6vJkMNsP+BLwFmnTx30jYrikDSRNLbQ7D7iZNO37HGn1z3m5bhHgStLU8Kukewe3joiPG/IJzMzMzHqg7jBNTES8S9o/sLz8ftKikdLrAA7Lj/K29wAr1C9KMzMzs96nWySDZgAjBq1YtW7FnRoYiJmZWRPpLtPEZmZmZtYFnAyamZmZNTEng2ZmZmZNzMmgmZmZWRPzAhIzMzOr2dAnx0Y9+9/560v4sIgG88igmZmZ9UiSRknarKxsY0ljy8q+K+leSVMkvSPpKUmHS5on17dI+ljS1MJjYuH6kDQtl78h6TRJszfkQzaAk0EzMzPrtST9GLgGuAJYOiIWAnYElgCWLDS9MiL6Fh79y7paPSL6Ahvl639a/+gbw9PEZmZm1itJEnAacFxEXFAqj4gXgV/NSp8R8YqkB4E1OiXIbsDJoJmZWTc09MmxVet2/voSDYykR1uBNAJ4bWd1KGkQsAFwcmf12dU8TWxmZma91cL5z/GlAkl/lzRR0vuSdi+03SGXlx73lvX1H0nTgBHAMODsukbeQB4ZNDMz62HuvmfZqnWbbjKygZF0e+/kPxcFXgOIiJ0AJD0AFBeBXBURu7XS15rASODHwIlAH+DDzg64K3hk0MzMzHqrF4E3gO07o7NIrgIeBn7XGX12Bx4ZNDPrRCMGrVi1bsWdGhiIWfOYs7RFTPZpbhMRn0g6GLhA0mTSquKJwNeARTrwnicC/5Z0YkSMb7N1N+dk0MzMzGrWDTeF/kfZ6weLLyLiSkmTgN8AfyZN7b4OnA9cXWi6o6Ttyvr6akS8Vf6GEfGspPuAQ4GDOxZ+13MyaGZmZj1SRAyssd0/gX+2Ut8CtLRS/4UEOCK2qOW9ewLfM2hmZmbWxJwMmpmZmTUxJ4NmZmZmTcz3DJqZWY/T6qrtF0Y0MBKzns/JoJmZWZ2MPeL+VuuXOHGDBkViVp2nic3MzMyamJNBMzMzsybmaWIzMzNg1YtXrVr37J7PNjASs8ZyMmhmZtZFTt1x66p1ix1xbgMjqd3d9ywb9ex/001G1v2EE0kBLBcRr0g6F3gjIn5fp/faALgwIlaoR/+dwdPEZmZm1rQiYp96JYK5//uLiaCkUZLektSnULaXpGGF15HbzFEomzOXRaFsmKQPJE0tPNZtb4xOBs3MzMwaa3bgwDbavAcUj7zbIpeV2z8i+hYeD7c3GCeDZmZm1iNJWjGPjk2UNFzSD3L5EElnSbpV0hRJj0hatkofQyQdn59vLGmspIPzKNybkn5SaDu3pFMkvS7pv5LOlfSlNmLcWNLYsuI/AYdI6t/KpZcCexRe7wFc0tp7zSrfM2hmZk1j4BG3Vq2br/o+1q06a597qtZt23/OWevU2iRpTuBm4G/A5sD6wI2S1spNdiKNpv0HuBj4Qy5ry1eA+YHFge8C10i6ISLeA04ElgXWAD4GrgB+B/ymneE/DgwDDgGOrtLmBuBXOWEUsAHQAhzfzvdqk5NBMzPrXVrmb6XyioaFYXX3LaAvcGJEfALcI+kWYOdcf31EPAog6XLgtBr7/Rg4LiJmAP+QNBVYQdIjwN7AahHxbu73BNIPVXuTQUhJ5IOS/q9K/QekZHdHUjJ4Uy4rd4akU/LzVyNizfYG4mTQzMzMeqLFgDE5ESwZTRrRAxhfKH+flDjW4p2cCJZfOwCYF3hC+nTBs0j3/7VbRDyXk9cjgGpnKF4C/DG/z+FV2hwQERfOSgwlTgbNzMza0NLSUrVuABs2LhArGgcsKWm2QkK4FPASMLAO7/c2MB1YOSLe6KQ+jyFNY59apf5+YFEggAdIU9SdzgtIzMzMrCd6hDRqd1jedmVjYBvg7/V4s5xwXgD8WdKXASQtLul7HejzFeBK4IAq9UH6TD/Iz+vCI4NmZmZWs0ZsCl2LiPhI0jbA2aR79t4A9oiIFwrTuJ3tcNK9fv+WtHB+z3OA2zvQ53HA7tUqI2J4B/quiZNBsyYyYlD15ZIr1rLGzsysG8mJ0kYVygeXvR4GLFF4rUpty9vlsoGF5x8AR+ZHrTGWv/fAsvoxwDxlZRWz2TySWIx941rjaI2nic3MzMyamJNBMzMzsw6QdGTZkXClx21dHVstPE1sZmZm1gERcQJwQlfHMas8MmhmZmbWxJwMmpmZmTUxJ4NmZmZmTczJoJmZmVkTczJoZmZm1sScDJqZmZk1MSeDZmZmZk3MyaCZmZlZE3MyaGZmZtbEnAyamZmZNTEng2ZmZmZNzMmgmZmZWRNzMmhmZmbWxJwMmpmZmTUxJ4NmZmZmTczJoJmZmVkTczJoZmZm1sS6RTIoaUFJ10uaJmm0pF2qtJOkkyS9kx8nSVKhfg1JT0h6P/+5RsM+hJmZmVkP1C2SQeAs4CNgEWBX4BxJK1dotzewHbA6sBqwDfALAElzATcClwELABcDN+ZyMzMzM6ugy5NBSX2AHwG/jYipEfEAcBOwe4XmewKnRsTYiHgDOBUYnOs2BuYATo+IDyPiDEDAJnX+CGZmZmY9liKiawOQvg48GBHzFsoOATaKiG3K2k4CNo+IR/LrtYB7I2I+Sf+b67YotL8l159a4X33Jo00AqwAvNjJH62eFgbe7uogejl/x43h77n+/B3XX0/8jpeOiAFdHYR1D3N0dQBAX2ByWdkkYL4qbSeVteub7xssr2utHyLifOD8WQm4q0l6PCLW6uo4ejN/x43h77n+/B3Xn79j6+m6fJoYmAr0KyvrB0ypoW0/YGqk4c329GNmZmZmdI9k8CVgDknLFcpWB4ZXaDs811VqNxxYrbi6mLTIpFI/ZmZmZkY3SAYjYhpwHXCcpD6S1gO2BS6t0PwS4CBJi0taDDgYGJLrhgEzgQMkzS1p/1x+Tz3j7yI9cnq7h/F33Bj+nuvP33H9+Tu2Hq3LF5BA2mcQ+BvwXeAd4IiIuELSBsBtEdE3txNwErBXvvRC4PA8TVxajHIhsBIwAvhZRDzZ0A9jZmZm1oN0i2TQzMzMzLpGl08Tm5mZmVnXcTJoZp1G0ihJm0k6UtKFXR2PWa1KP7tdHYdZV3AyaL1W/uU+XdJUSf+VNERSX0nDJH2QyydJuk/SqoXrWiSFpAPL+jswl7cUyvpJOl3S67m/kfn1wg38qN1ORJwQEXu13dIqyT+7b+UTmkple0kalp9vK+kpSZMlvS3pHknLFNouL+nqXDdJ0jOSDpI0u6SB+ee44j6zrV1b9w/eBCQNlvRAV8dhVuRksBfxL+uKtskLkNYE1gKOzuX75/IFSSvRy1evvwTsUVa2Zy4HPj0P+25gZeD7pH0t1yUtgvpmp34KA5ruZ3x24MDyQklfI+2scDAwP7AM6Xz3mbl+WeARYAywakTMD/yY9PNfcRP+Qt+zfK21rVoCbtbVnAx2ovyv+d9Iel7Se5IukjSPpAUk3SJpQi6/RdISheuGSfq9pAclTZF0R3FkSdL6kh6SNFHSGEmDc/kQSedI+oekacB3JC0m6dr8Xq9JOqDQzzclPZz7eVPSmTmhQcmf82jEZEnPSlqlcd9efeWzrG8DVikrnwn8nbQCvegxYF5JKwPkP+fJ5SV7AEsBP4yI5yPik4h4KyJ+HxH/qNNH6RHy6Opl+flt+myrp1L905K2z88HSbpT0ruSXpS0Q6HdF37GG/pButafgEMk9S8rXwN4LSLujmRKRFwbEa/n+mOBhyLioIh4EyAiXoyIXSJiYhvv2ZFre4M18kjoJElXSpoHQNLWeSR2Yv5dvFrpAklH5BmBKfl3/w8LdYPz7/U/S3oHuBI4F1hXaSZhYqM/oFklTgY7367A94BlgeVJI1GzARcBS5OSh+nAmWXX7QL8BPgyMBdwCICkpUlJzF+AAaS/CJ4qu+4PpH+1PwTcDDwNLA5sCvxa0vdy25nA/5LO0Vw31++X6zYHNswxzw/sQBrh6hUkLQlsCTxZVj4X6b/ZvytcdimfjQ7uyRdHDzcD/hkRUzs32l5nKLBz6YWklUj/L9yqNA16J3AF6Wd/J+Ds3Kak+DPeTNNrj5NGrQ8pK/8PMCgnGN+R1LesfjPgmll8z45c2xvsQBrlX4Z0aMFgpS3L/gb8AlgIOA+4SdLc+ZqRwAak35vHApdJWrTQ5zrAq8AiwG7APsDDEdE3IvrX/ROZ1cDJYOc7MyLGRMS7pL/Ado6Id/K/3N+PiCm5fKOy6y6KiJciYjpwFSnpg/QX4V0RMTQiPs59PVW47saIeDAiPgFWBQZExHER8VFEvApcQPoLloh4IiL+HREzImIU6ZdaKY6PSX/ZDiJtOTSiNDLQw92Q//X9APAv4IRcfkYunwLsT/olXu4yYGdJc5K+w8vK6hcCesN3VG/Xk0Zcls6vdwWui4gPga2BURFxUf65fBK4ljQ1WfLpz3hEfNDY0Lvc74BfSRpQKsj/X29M+gffVcDbeQS1lBR25Oey2X+mz4iIcfn3982k38N7A+dFxCMRMTMiLgY+BL4FEBFX52s+iYgrgZf5/G0i4yLiL/nne3pjP45ZbZwMdr4xheejgcUkzSvpPEmjJU0G7gP6l93/NL7w/H2g9It9SdK/PGt5v6Xz+00sPYAjSf8iLd0Yfouk8TmOE0ijhETEPaTRyrOAtySdL6n8rOeeaLuI6B8RS0fEfoVfxgfkf5V/iZSQXFOc+gHI026vkL6nlyOi+F1DGjldFGtV/gfQreR/lJBGCS/Pz5cG1in7md0V+Eqhi/LvvWlExHPALcARZeX/jogdImIAaVRqQ+CoXN2Rn8tm/5mu9Ht4aeDgsp/RJYHFACTtUZhCnki6FaW4gKxpf36t53Ay2PmWLDxfChhHutF7BWCdiOhH+sUNINo2hjTlXE1x1/AxpHuJ+hce80XElrn+HOAFYLkcx5HFGCLijIj4Bun+ueWBQ2uIr0fL/5q/n5T0bV6hSelG/Usq1N0FfE+FFZ9W1VDSKOu6pHsv783lY4B/lf3M9o2IfQvXNvvO+McAPyeNBH5BRDxGOtKzdD/sXcCPZvG9OnJtbzUG+EPZz+i8ETE0j3ZfQJpdWCj/A/M5Pv+7vfznt9l/nq0bcjLY+X4paQmlI/aOIt0wPB/pPsGJufyYdvR3ObCZpB0kzSFpIUlrVGn7KDBF0uGSvqS0jcQqktbO9fMBk4GpkgYBn/6FK2ltSevkKdFpwAfAJ+2Is8fKCcpKwPAK1VeSksSrKtRdSvqL4tq8CGK2/N/nSElbVmjfzP5BGmE5Drgy39YAadRreUm7S5ozP9aWtGKXRdrNRMQrpJ/DA+DTBWU/l/Tl/HoQ8AM+u+/1GODbkv4k6Su5zdckXVa2GGVupQVupcds7bi2mVwA7JN/P0pSH0lbSZoP6ENK7iYASPoJZYvUKvgvsES+X9msW3Ay2PmuAO4g3TA8EjgeOJ00Hfk26Rf2P2vtLE9VbkkanXqXtHhk9SptZ5KmPNcAXsvvdyHpxmZIN6LvQrpP7gLSXzAl/XLZe6Tp7XdIqxl7qzPzar6ppKTu6Ii4rbxRREyPiLsq3euT73nbjDTaeicp0X6UNEX0SF2j72Hyd3Ud6fu6olA+hZRs70QaRR9POn987grdNLPjSIkHwERS8vds/vn9J+m+zJMBImIkaYHYQGC4pEmk+zAfJ/2/XzKV9I/U0mOTdlzbNCLicdLI7Jmk34+vAINz3fPAqcDDpCRvVeDBNrq8h/QPz/GS3q5P1Gbt47OJO5GkUcBeEXFXV8diZmZmVguPDJqZmZk1MSeDZmZmZk3M08RmZmZmTcwjg2ZmZmZNzMmgmZmZWRNzMmhmZmbWxJwMmpmZmTUxJ4NmZmZmTczJoJmZmVkTczJoZmZm1sScDJqZmZk1MSeDZmZmZk3MyaCZmZlZE3MyaGZIGiIp8uNjSW9JulfSLyXN2Y5+Ns59LFzPeCu878D8vms18n3NzHoDJ4NmVnIXsCgwENgcuBk4FrhfUp8ujMvMzOrIyaCZlXwYEeMj4o2IeCoiTgM2BtYEDgOQtJukxyRNyaOHV0taPNcNBO7NfU3II3VDct33Jd0v6T1J70q6XdKKxTeX9DtJoyV9KGm8pEsKdZJ0mKSRkqZLelbSboXLX8t/Ppbfd1i+blVJd0uaLGmqpKclfaeTvzczsx7NyaCZVRURzwH/BH6Ui+YCjgFWB7YGFgaG5roxhXYrk0YZD8yv+wCnA98kJZiTgJslzQUg6UfAIcB+wHK570cLoRwP/Az4JbAS8EfgPElb5fpv5j+/n993+/z6CuDNXL8G0AJ80N7vwcysN5ujqwMws27veWAzgIj4W6H8VUn7AiMkLRERYyW9m+veioi3Sw0j4tpih5J+AkwmJWkPAEuTkrY7IuJj4HXg8dy2D3AQsHlE3J+7eE3SN0nJ4a3AhFz+TkSML7zV0sApEfFCfv3KrH4JZma9lUcGzawtAgJA0pqSbszTuVPICRuwVKsdSMtKuiJP804G/kv6/VO67mpgHlKS91dJP5Y0d65bKdf9M0/1TpU0FdgXWLaN2E8DLpR0j6SjJA1q30c3M+v9nAyaWVtWIo0C9gFuB94HdgfWJk3LQpo+bs0twADgF8A6wNeBGaXrImIMsEKunwycCjyR37P0e2ob0lRv6bEyaaFLVRHRkuO/Afg28Iykn7YRq5lZU/E0sZlVJWkVUsJ3PDCIdI/gkRHxWq7fvuySj/Kfsxf6WChfu19E3JvL1qTs909EfECa8r1V0onAeGA94GHgQ2DpiLinSqhfeN9Cvy8DLwNnSDoH2Av4W3k7M7Nm5WTQzErmlvQV0kjcAGBT4EjgCeAUYF5SUra/pLOAFYHfl/UxmjSlvJWkm4HpwHvA28DPJY0BFgf+RBoZBEDSYNLvo0eAqcCOwMfAyxExRdIpwCmSBNwH9AW+BXwSEecDb+X3+p6kUaRFIh/luK8GRgGLAOvn9zAzs8zTxGZWshlpEcfrwN3AD0irbzeMiGkRMQHYE9iOtKjkGNLCjk9FxBu5/A+k+wLPjIhPSMndasBzwFnAb0mJZclE0mrh+3ObHwHbl0Ygc/sW0orj4cCduc1r+X1nAAeQRv3GATcCM4EFgCHAi8D1pFHGz8VsZtbsFBFdHYOZmZmZdRGPDJqZmZk1MSeDZmZmZk3MyaCZmZlZE3MyaGZmZtbEnAyamZmZNTEng2ZmZmZNzMmgmZmZWRNzMmhmZmbWxP4fvPnW0QotwvAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# score = pd.read_csv('~/SCALEX/notebook/benchmark/overcorrect_benchmark.txt', sep='\\t', index_col=0)\n",
    "\n",
    "datasets = list(score.index)\n",
    "methods = list(score.columns)\n",
    "palette = { 'SCALEX': (0.8392156862745098, 0.15294117647058825, 0.1568627450980392),\n",
    "            'Raw': (0.09019607843137255, 0.7450980392156863, 0.8117647058823529),\n",
    "            'Seurat_v3': (1.0, 0.4980392156862745, 0.054901960784313725),\n",
    "            'Harmony':  (0.12156862745098039, 0.4666666666666667, 0.7058823529411765),\n",
    "            'Conos':(0.4980392156862745, 0.4980392156862745, 0.4980392156862745),\n",
    "            'BBKNN': (0.5803921568627451, 0.403921568627451, 0.7411764705882353),\n",
    "            'Scanorama': (0.8901960784313725, 0.4666666666666667, 0.7607843137254902),\n",
    "            'FastMNN':  (0.17254901960784313, 0.6274509803921569, 0.17254901960784313),\n",
    "            'scVI': (0.5490196078431373, 0.33725490196078434, 0.29411764705882354),\n",
    "            'online_iNMF':(0.7372549019607844, 0.7411764705882353, 0.13333333333333333),\n",
    "            'LIGER':(0.6509803921568628, 0.807843137254902, 0.8901960784313725),}\n",
    "\n",
    "score = pd.DataFrame({'Dataset':np.repeat(datasets, len(methods)),\n",
    "                      'Method': methods* len(datasets),\n",
    "                      'overcorrect_score': np.reshape(score.values, (1, -1)).squeeze(),})\n",
    "score = score[score['Method'] !='Raw']\n",
    "legend_params = {'loc': 'center left',\n",
    "                 'bbox_to_anchor':(1.01, 0.35),\n",
    "                 'fontsize': 12,\n",
    "                 'ncol': 1,\n",
    "                 'frameon': False,\n",
    "                 'markerscale': 1\n",
    "                }\n",
    "\n",
    "plt.rcParams['savefig.dpi'] = 100\n",
    "plt.rcParams['figure.dpi'] = 100\n",
    "plt.figure(figsize=(8,6))\n",
    "fig = sns.barplot(x='Dataset', y='overcorrect_score', data=score, hue='Method',palette=palette, saturation=1)\n",
    "plt.ylabel('Over-correctiont Score',fontsize=14, labelpad=10)\n",
    "plt.xlabel('Datasets',fontsize=14, labelpad=10)\n",
    "plt.xticks(fontsize=12)\n",
    "plt.yticks(fontsize=12)\n",
    "plt.title('',fontsize=16,y=1.02)\n",
    "plt.legend(**legend_params)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
